rm(list = ls())
gc()
options(stringsAsFactors = F)
library(tidyverse)
library(clusterProfiler)
library(msigdbr)  #install.packages("msigdbr")
library(GSVA) 
library(GSEABase)
library(pheatmap)
library(limma)
library(BiocParallel)
try({
  setwd('../gsva/')
})

try({
  setwd('./gsva/')
})

## msigdbr包提取下载 一般尝试KEGG和GO做GSVA分析
##KEGG
KEGG_df_all <-  msigdbr(species = "Homo sapiens", # Homo sapiens or Mus musculus
                        category = "C2",
                        subcategory = "CP:KEGG") 
KEGG_df <- dplyr::select(KEGG_df_all,gs_name,gs_exact_source,gene_symbol)
kegg_list <- split(KEGG_df$gene_symbol, KEGG_df$gs_exact_source) ##按照gs_name给gene_symbol分组

##GO
GO_df_all <- msigdbr(species = "Homo sapiens",
                     category = "C5")  
GO_df <- dplyr::select(GO_df_all, gs_name, gene_symbol, gs_exact_source, gs_subcat)
GO_df <- GO_df[GO_df$gs_subcat!="HPO",]

# 筛选 "gs_subcat" 列中第4-5个字符为 "BP" 的行
BP_df <- GO_df[substr(GO_df$gs_subcat, 4, 5) == "BP", ]
# 筛选 "gs_subcat" 列中第4-5个字符为 "BP" 的行
CC_df <- GO_df[substr(GO_df$gs_subcat, 4, 5) == "CC", ]
# 筛选 "gs_subcat" 列中第4-5个字符为 "BP" 的行
MF_df <- GO_df[substr(GO_df$gs_subcat, 4, 5) == "MF", ]
go_list <- split(BP_df$gene_symbol, BP_df$gs_exact_source) ##按照gs_name给gene_symbol分组
##提取自己的样本
FPKM=read.table("brca_counts_MRNA.txt",sep="\t",header=T,check.names=F)
dup_rows <- FPKM[duplicated(FPKM[, 1]), 1]
dup_suffix <- ave(dup_rows, dup_rows, FUN = seq_along)
FPKM[duplicated(FPKM[, 1]), 1] <- paste0(dup_rows, ".", dup_suffix)
rownames(FPKM) <- FPKM[,1]
FPKM <- FPKM[,-1]
write.csv(FPKM,"BRCACOUNTS自测.csv", sep = ",",row.names = T,col.names = NA,quote = F)

####  GSVA  ####
#GSVA算法需要处理logCPM, logRPKM,logTPM数据或counts数据的矩阵####
#dat <- as.matrix(counts)
dat <- as.matrix(log2(edgeR::cpm(FPKM))+1)
#dat <- as.matrix(log2(tpm+1))

#dat <- as.matrix(log2(FPKM+1))
##这一步是提取自己的基因集
#genelist <- read.table("工作簿3.csv",sep=",",header=T,check.names=F)
#geneset <- lapply(genelist, function(col) col[!is.na(col) & col != ""])

geneset <- go_list##这一步是用go的
geneset <- kegg_list
##linux系统

huangPar2 <- gsvaParam(exprData = dat, geneSets = geneset,kcdf = "Gaussian")
gsva_mat2 <- gsva(huangPar2)#调用所有核


gsva_mat <- gsva(expr=dat, 
                 gset.idx.list=gene.set, 
                 kcdf="Gaussian" ,#"Gaussian" for logCPM,logRPKM,logTPM, "Poisson" for counts
                 verbose=T, 
                 parallel.sz = parallel::detectCores())#调用所有核
##win系统

#param <- SnowParam(workers = parallel::detectCores(), type = "SOCK")
##gsva_mat <- gsva(expr = dat, gset.idx.list = geneset, kcdf = "Gaussian", BPPARAM = param)
##CPU会直接占满卡死，不能用
counts=read.table("CGGA_FPKM.txt",sep="\t",header=T,check.names=F, row.names = 1)
####  GSVA  ####
#GSVA算法需要处理logCPM, logRPKM,logTPM数据或counts数据的矩阵####
#dat <- as.matrix(counts)
#dat <- as.matrix(log2(edgeR::cpm(counts))+1)
#dat <- as.matrix(log2(tpm+1))
dat <- as.matrix(log2(counts+1))
##这一步是提取自己的基因集
#genelist <- read.table("工作簿3.csv",sep=",",header=T,check.names=F)
#geneset <- lapply(genelist, function(col) col[!is.na(col) & col != ""])

#geneset <- go_list##这一步是用go的

##linux系统

huangPar2 <- gsvaParam(exprData = dat, geneSets = geneset,kcdf = "Gaussian",minSize =1 )

gsva_mat2 <- gsva(expr=dat, 
                 gset.idx.list=geneset, 
                 kcdf="Gaussian" ,#"Gaussian" for logCPM,logRPKM,logTPM, "Poisson" for counts
                 verbose=T, 
                 parallel.sz = parallel::detectCores())#调用所有核

gsva_mat2 <- gsva(huangPar2,verbose = TRUE)
##win系统
write.csv(gsva_mat,"tcga_brca_kegg_matrix.csv", sep = ",",row.names = T,col.names = NA,quote = F)
write.csv(gsva_mat2,"CGGA_gsva_go_matrix.csv", sep = ",",row.names = T,col.names = NA,quote = F)

#### 进行limma差异处理 ####
##设定 实验组exp / 对照组ctr
gl
library(limma)
group_list <- read.table("Untitled",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
group_list <- as.matrix(group_list)
exp="low"
ctr="high"
gsva_mat <- read.table("TCGA_gsva_go_matrix.csv",sep = ",",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(gsva_mat)
contrast.matrix <- makeContrasts(contrasts=paste0(exp,'-',ctr),  #"exp/ctrl"
                                 levels = design)

fit1 <- lmFit(gsva_mat,design)                 #拟合模型
fit2 <- contrasts.fit(fit1, contrast.matrix) #统计检验
efit <- eBayes(fit2)                         #修正

summary(decideTests(efit,lfc=1, p.value=0.05)) #统计查看差异结果
tempOutput <- topTable(efit, coef=paste0(exp,'-',ctr), n=Inf)
degs <- na.omit(tempOutput) 
write.csv(degs,"gsva_go_degs.results.csv")



